﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace fftwlib
{
    public class OctaveCalculate
    {
        /// <summary>
        /// 通过功率谱计算1/3倍频程
        /// </summary>
        /// <param name="powerSpectrum">功率谱</param>
        /// <param name="freqRes">频率分辨率</param>
        /// <param name="lowdB">低频段总值</param>
        /// <param name="highdB">高频段总值</param>
        /// <param name="totaldB">全频段总值</param>
        /// <returns>1/3倍频程</returns>
        public double[] OneThirdOctave_10Hz_10000Hz(float[] powerSpectrum, float freqRes, ref double lowdB, ref double highdB, ref double totaldB)
        {
            double[] tempOctave = new double[31];
            double[] _1_3octavedB = new double[31];
            double c = Math.Pow(2, 1.0 / 3.0);
            for (int m = 0; m < 31; m++)
            {//开始计算1/3倍频程
                double upLeft, downRight = 0;
                double fl = 10 * Math.Pow(c, m - 0.5);
                double fu = 10 * Math.Pow(c, m + 0.5);
                int n = (int)Math.Ceiling(fl / freqRes);
                int t = (int)Math.Ceiling(fu / freqRes) + 1;
                downRight = (freqRes * (float)n + freqRes * (float)(n + 1)) / 2.0;
                upLeft = (freqRes * (float)t + freqRes * (float)(t - 1)) / 2.0;
                if (downRight < fl)
                {
                    n++;
                    downRight = (freqRes * (float)n + freqRes * (float)(n + 1)) / 2.0;
                }
                if (upLeft > fu)
                {
                    t--;
                    upLeft = (freqRes * (float)t + freqRes * (float)(t - 1)) / 2.0;
                }

                if (n == t)
                {
                    tempOctave[m] += Math.Pow((fu - fl) / freqRes * powerSpectrum[n], 2);
                }
                else if (t == n + 1)
                {
                    tempOctave[m] += Math.Pow((downRight - fl) / freqRes * powerSpectrum[n], 2) + Math.Pow((fu - upLeft) / freqRes * powerSpectrum[t], 2);
                }
                else
                {
                    tempOctave[m] += (Math.Pow((downRight - fl) / freqRes * powerSpectrum[n], 2) + Math.Pow((fu - upLeft) / freqRes * powerSpectrum[t], 2));
                    for (int o = n + 1; o < t; o++)
                    {
                        tempOctave[m] += Math.Pow(powerSpectrum[o], 2);
                    }
                }
            }

            //将1/3倍频程单位从m/s2转换为dB,同时计算振动总值
            for (int n = 0; n < 31; n++)
            {
                _1_3octavedB[n] = 20.0 * (6.0 + Math.Log10(Math.Sqrt(Math.Abs(tempOctave[n]))));

                if (n < 16)
                {//计算低频段总值
                    lowdB += Math.Pow(10.0, _1_3octavedB[n] / 10.0);
                }
                else
                {//高频段
                    highdB += Math.Pow(10.0, _1_3octavedB[n] / 10.0);
                }

                totaldB += Math.Pow(10.0, _1_3octavedB[n] / 10.0);
            }
            lowdB = 10.0 * Math.Log10(lowdB);
            highdB = 10.0 * Math.Log10(highdB);
            totaldB = 10.0 * Math.Log10(totaldB);

            return _1_3octavedB;
        }
    }
}
